Prostate Cancer ~ Predictdb eQTL&sQTL from Prostate & Testis & Lymphocytes
Nina Yifan Yang
2024-10-15
Last updated: 2024-10-15
Checks: 7 0
Knit directory: Multigroup_cTWAS/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20241004) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 741912e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: analysis/.Rhistory
Ignored: analysis/figure/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown
(analysis/cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd)
and HTML
(docs/cTWAS_noLD_european_prostatetestislymphocytes_output.html)
files. If you’ve configured a remote Git repository (see
?wflow_git_remote), click on the hyperlinks in the table
below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 741912e | yy896 | 2024-10-15 | wflow_publish("analysis/cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | 3355d2f | yy896 | 2024-10-07 | Build site. |
| Rmd | 01306bd | yy896 | 2024-10-07 | wflow_publish("cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | 735a51d | yy896 | 2024-10-07 | Build site. |
| Rmd | fb77724 | yy896 | 2024-10-07 | wflow_publish("cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | 9d60379 | yy896 | 2024-10-07 | Build site. |
| Rmd | e71cea7 | yy896 | 2024-10-07 | wflow_publish("cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | 314c305 | yy896 | 2024-10-07 | Build site. |
| Rmd | cfdba1a | yy896 | 2024-10-07 | wflow_publish("analysis/cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | 93f8704 | yy896 | 2024-10-07 | Build site. |
| Rmd | a480298 | yy896 | 2024-10-07 | wflow_publish("analysis/cTWAS_noLD_european_prostatetestislymphocytes_output.Rmd") |
| html | fce136d | yy896 | 2024-10-04 | Build site. |
| Rmd | a715301 | yy896 | 2024-10-04 | Add my first analysis |
We present a sample cTWAS report based on real data analysis. The analyzed trait is prostate cancer, the prediction models are Prostate & Testis & Lymphocytes gene expression and splicing models trained on GTEx v8 in the PredictDB format.
Analysis Settings
Input Data
- GWAS Z-scores
The summary statistics for prostate cancer are downloaded from https://www.ebi.ac.uk/gwas/downloads/summary-statistics, using dataset ID: GCST90274714.
The sample size is
[1] "gwas_n = 726828"- Prediction models
The prediction models used in this analysis are prostate & testis & lymphocytes gene expression and splicing models, trained on GTEx v8 in the PredictDB format. These models can be downloaded from https://predictdb.org/post/2021/07/21/gtex-v8-models-on-eqtl-and-sqtl/
- Reference data
We did not include LD reference in the analysis.
Data processing and harmonization
We map the reference SNPs to regions following the instructions from the cTWAS tutorial.
When processing z-scores, we exclude multi-allelic and
strand-ambiguous variants by setting
drop_multiallelic = TRUE and
drop_strand_ambig = TRUE.
The process can be divided into steps below, users can expand the code snippets below to view the exact code used.
- Input and output settings
weight_prostate_expression_file <- system.file("extdata/project_data",
"mashr_Prostate_e.db", package = "ctwas")
weight_prostate_splicing_file <- system.file("extdata/project_data",
"mashr_Prostate_s.db", package = "ctwas")
weight_testis_expression_file <- system.file("extdata/project_data",
"mashr_Testis_e.db", package = "ctwas")
weight_testis_splicing_file <- system.file("extdata/project_data",
"mashr_Testis_s.db", package = "ctwas")
weight_lymphocytes_expression_file <- system.file("extdata/project_data",
"mashr_Cells_EBV-transformed_lymphocytes_e.db", package = "ctwas")
weight_lymphocytes_splicing_file <- system.file("extdata/project_data",
"mashr_Cells_EBV-transformed_lymphocytes_s.db", package = "ctwas")
z_snp <- vroom::vroom("~/Desktop/Prostate Cancer Project/GCST90274714.mod.tsv.gz", col_names = TRUE, show_col_types = FALSE)
z_snp$Z <- z_snp$beta / z_snp$standard_error
z_snp <- z_snp[, c("rsid", "other_allele", "effect_allele", "Z")]
colnames(z_snp) <- c("id", "A1", "A2", "z")
region_file <- system.file("extdata/ldetect", "EUR.b38.ldetect.regions.RDS", package = "ctwas")
region_info <- readRDS(region_file)
## other parameters
ncore <- 6- Preprocessing GWAS
### Preprocess SNP_map
ref_snp_info_file <- system.file("extdata/project_data", "ukb_b38_0.1_var_info.Rvar.gz", package = "ctwas")
ref_snp_info <- data.table::fread(ref_snp_info_file, sep = "\t")
class(ref_snp_info) <- "data.frame"
res <- create_snp_map(region_info, ref_snp_info)
region_info <- res$region_info
snp_map <- res$snp_map
### Preprocess GWAS z-scores
z_snp <- preprocess_z_snp(z_snp, snp_map,
drop_multiallelic = TRUE,
drop_strand_ambig = TRUE)- Preprocessing weights
weights_prostate_expression <- preprocess_weights(weight_prostate_expression_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "expression",
context = "prostate",
weight_name = "prostate_expression",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights_prostate_splicing <- preprocess_weights(weight_prostate_splicing_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "splicing",
context = "prostate",
weight_name = "prostate_splicing",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights_testis_expression <- preprocess_weights(weight_testis_expression_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "expression",
context = "testis",
weight_name = "testis_expression",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights_testis_splicing <- preprocess_weights(weight_testis_splicing_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "splicing",
context = "testis",
weight_name = "testis_splicing",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights_lymphocytes_expression <- preprocess_weights(weight_lymphocytes_expression_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "expression",
context = "lymphocytes",
weight_name = "lymphocytes_expression",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights_lymphocytes_splicing <- preprocess_weights(weight_lymphocytes_splicing_file,
region_info,
gwas_snp_ids = z_snp$id,
snp_map = snp_map,
type = "splicing",
context = "lymphocytes",
weight_name = "lymphocytes_splicing",
weight_format = "PredictDB",
drop_strand_ambig = TRUE,
scale_predictdb_weights = TRUE,
load_predictdb_LD = TRUE,
filter_protein_coding_genes = TRUE,
ncore = 6)
weights <- c(
weights_prostate_expression,
weights_prostate_splicing,
weights_testis_expression,
weights_testis_splicing,
weights_lymphocytes_expression,
weights_lymphocytes_splicing)Running cTWAS analysis
We use the ctwas main function ctwas_sumstats_noLD to
run the cTWAS analysis without LD. For more details on this function,
refer to the cTWAS tutorial: https://xinhe-lab.github.io/multigroup_ctwas/articles/minimal_tutorial.html
All arguments are set to their default values.
Users can expand the code snippets below to view the exact code used.
ctwas_res <- ctwas_sumstats_noLD(z_snp,
weights,
region_info,
snp_map,
thin = 0.1,
ncore = 6)Parameter estimation
ctwas_res is the object contains the outputs of
cTWAS
We extract the estimated parameters by
param <- ctwas_res$param
We make plots using the function
make_convergence_plots(param, gwas_n) to see how estimated
parameters converge during the execution of the program:
| Version | Author | Date |
|---|---|---|
| fce136d | yy896 | 2024-10-04 |
These plots show the estimated prior inclusion probability, prior effect size variance, enrichment and proportion of variance explained (PVE) over the iterations of parameter estimation. The enrichment is defined as the ratio of the prior inclusion probability of molecular traits over the prior inclusion probability of variants. We generally expect molecular traits to have higher prior inclusion probability than variants. Enrichment values typically range from 20 - 100 for expression traits.
Then, we use summarize_param(param, gwas_n) to obtain
estimated parameters (from the last iteration) and to compute the PVE by
variants and molecular traits.
[1] "The number of genes/introns/SNPs used in the analysis is:" prostate|expression prostate|splicing testis|expression
9693 29875 11361
testis|splicing lymphocytes|expression lymphocytes|splicing
50036 8478 26609
SNP
7146350 ctwas_parameters$prop_heritability contains the
proportion of heritability mediated by molecular traits and variants, we
visualize it using pie chart.
ctwas_parameters <- summarize_param(param, gwas_n)
data <- data.frame(
category = names(ctwas_parameters$prop_heritability),
percentage = ctwas_parameters$prop_heritability
)
# Calculate percentage labels for the chart
data$percentage_label <- paste0(round(data$percentage * 100, 1), "%")
ggplot(data, aes(x = "", y = percentage, fill = category)) +
geom_bar(stat = "identity", width = 1) +
coord_polar("y", start = 0) +
theme_void() + # Remove background and axes
geom_text(aes(label = percentage_label),
position = position_stack(vjust = 0.5)) +
scale_fill_brewer(palette = "Pastel1") +
labs(fill = "Category") +
ggtitle("Percent of heritability")
| Version | Author | Date |
|---|---|---|
| fce136d | yy896 | 2024-10-04 |
Diagnosis plots
For all genes analyzed, we compare the z-scores and fine-mapping PIPs”
ggplot(data = finemap_res[finemap_res$type!="SNP",], aes(x = abs(z), y = susie_pip)) +
geom_point() +
labs(x = "abs(z-scores)", y = "PIPs") +
theme_minimal()
| Version | Author | Date |
|---|---|---|
| fce136d | yy896 | 2024-10-04 |
Fine-mapping results
We process the fine-mapping results here.
We first add gene annotations to cTWAS results.
mapping_table is the combination of both expression and
splicing data.
finemap_res <- anno_finemap_res(finemap_res,
snp_map = snp_map,
mapping_table = mapping_table,
add_gene_annot = TRUE,
map_by = "molecular_id",
drop_unmapped = TRUE,
add_position = TRUE,
use_gene_pos = "mid")The annotated fine-mapping results, ones within credible sets are shown
finemap_res_show <- finemap_res[!is.na(finemap_res$cs) &finemap_res$type !="SNP",]
DT::datatable(finemap_res_show,caption = htmltools::tags$caption( style = 'caption-side: topleft; text-align = left; color:black;','The annotated fine-mapping results, ones within credible sets are shown'),options = list(pageLength = 5) )Warning in instance$preRenderHook(instance): It seems your data is too big for
client-side DataTables. You may consider server-side processing:
https://rstudio.github.io/DT/server.html